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Applications  of  Differential  Topology  to  Grid  Generation 


D.C.  Wilson 
University  of  Florida 


Abstract 


)  The  purpose  of  this  paper  is  to  indicate  how  smoothing 
techniques  from  Differential  Topology  can  be  applied  to 
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the  area  of^grid  generation  in  Computational  Fluid  Dynamics 
The  basic  method  is  to  patch  together  one  global  grid  from 
a  number  of  smaller  ones.  The  smoothing  theory  allows  one 


to  blend  the  grid  from  one  section  into  the  grid  of  an 
adjacent  one.  — 
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I .  Introduction 


The  author  developed  the  ideas  outlined  in  this  paper  while 
attempting  to  construct  an  algebraic  grid^-  for  the  X-24C  aircraft. 
(See  Figure  1.)  In  a  natural  way  this  plane  can  be  divided 
into  the  forebody,  the  body,  and  the  airfoil.  The  forebody 
can  be  further  divided  into  the  nose  and  the  canopy.  The  airfoil 
region  can  be  subdivided  into  the  airfoil  and  the  body  adjacent 
to  it.  The  frustration  of  trying  to  blend  these  various  pieces 
together  into  one  large  grid  drove  the  author  to  search  for 
a  technique  of  approximating  one-to-one  continuous  functions 
(i.e.  homeomorphisms )  by  one-to-one  smooth  functions  (i.e.  diffeo- 
morphisms).  Techniques  lifted  from  the  Differential  Topology 
books  by  Hirsch2  and  Munkres^  provide  the  foundation  for  the 
theorems  and  ideas  presented  in  Section  II.  Theorems  2.1  and 
2.3  guarantee  that  any  continuous  transformation  can  be  approximated 
arbitrarily  closely  by  a  smooth  one.  Theorems  2.4  and  3.1 
can  be  used  to  ensure  that  the  approximation  will  have  non-zero 
Jacobian.  Theorem  2.5  ensures  that  the  approximation  of  an 
orthogonal  grid  will  be  almost  orthogonal. 

To  apply  the  smoothing  theory  to  grid  generation  a  large 
number  of  double  (or  triple)  integrals  must  be  computed.  Theorem 
3.2  is  the  2-dimensional  version  of  Simpson's  Rule  which  was 
used  to  calculate  the  convolutions  at  the  various  grid  points 
for  the  examples  pictured  at  the  end  of  the  paper.  While  other 

more  sophisticated  methods  (e.g.  iterated  Romberg)  could  be 
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II.  The  Smoothing  Theory 


An  excellent  reference  for  the  Advanced  Calculus  needed  in 
this  section  is  Rudin^.  While  the  theorems  are  all  stated  for 
two  dimensions,  they  are  all  valid  in  n  dimensions. 

The  smoothing  techniques  from  Differential  Topology  involve 
a  convolution  with  a  "bump"  function.  While  there  are  innumerable 
choices  for  a  bump  function  (e.g.  cubic  B-splines),  the  one 
indicated  below  proved  convenient. 


Definition . 


a  (x)  = 


0  if  x  <  0 

e-l./x  if  x  >  o 


The  real  valued  function  a  is  C“  (i.e.  All  higher  derivatives 
exists);  the  graph  of  a  is  indicated  in  Figure  2. 

If  a  <  b,  then  define  0(x)  -  Q!(x  -  a)  •  cx(b  -  x)  .  The 
function  B  is  also  C°°;  its  graph  is  indicated  in  Figure  3. 

Since  we  want  the  interval  [a,b]  to  be  symmetric  about  zero, 

EXP(-p2/(p2  -  )  if  |  x|  <  p 


we  define 


Bp(x)  = 


if  |  x|  >  Pi  where  the 


parameter  p  >  0.  If  A  =  f  B  (x)dx  and  0  (x)  =  r  (x)/A  .  then 

P  -p  P  P  P  P 

0p(x)  is  a  C  bump  function  such  that  0p(x)dx  =  1.  The 
function  0p(x,y)  =  0p(x)  •  Op(y)  is  a  bump  function  of  two 


variables  such  that  0  (x,y)dxdy  =  1.  I 

-p-p  P 

Let  Rp  denote  the  square  £-p,p|  x  £-P'P]  * 

Def inition .  If  F(^,n)  is  a  piecewise  continuous  function  from 
2 

R  -*■  R,  the  convolution  of  F  by  0  is: 


0  *  F(r,n)  =  /  J  F ( r -u  ,  n  -v  )  *Q  ( u ,  v  )dudv . 

P  4  1  Rp  4  p 

This  convolution  can  be  thought  of  as  an  "average"  of  the  values 
F  near  Heuri stica lly /  Theorem  2.1  states  that  if  p 

is  "small",  then  0p*  F  is  "close"  to  F. 

Theorem  2.1. 

If  F  is  continuous,  e  >0,  and|F(£-u,  n~v)  -  F(£,rj)l<  c 
for  all  (u,v)  £  R^,  then  |0p*F(^,n)  ~  F(£,n)|<  e  . 

Proof . 

The  result  follows  from  the  following  sequence  of 
equalities  and  inequalities  * 

I 0p*  F(£,n)  -  F(C»n)| 

=  I JJ  F(£-u,  n-v)0  (u,v)dudv  -  F ( c , n ) | 

Rp  p 

<  JJ  0  (u,v) • |F( C-u,  n-v)  -  F(5,n)j  dudv 
Rp 

5  eJJ  0  (u,v)dudv  =  t. 

Rp  P 

Theorem  2.2  shows  that  if  p  is  "small”  and  9F ( 5 , n )  is 

is  continuous  at  (  £  ,  n  )  then  i Q-P..*.  .F '  n  —  is  "close" 

9  £ 

to  9F ( 5 , n  )  .  Note  that  Theorem  2.2  can  be  easily  generalized  to 
95 

higher  derivatives. 

Theorem  2.2. 

If  F( C , n )  is  G1  ,  then  (0p  »  F( 5>n)  ) 


o  _  * 


9F( 5 , n ) ) 


FI-' 


«*_"  v* 


ipi  r  m 


5. 


v1 


Moreover,  if  c  >  0  and  1 3F ( C-u ,  n-v )  -  F ( C , n ) 1  -  c  for  (u,v)  £  Rp, 

H  at 


then  | 30p  *  F( t , n )  -  3F(t#n 


<  E 


3£ 


3t 


3® 


Proof . 


I 


k 


Since  jj  0p ( u , v )f( t~u ,  n-v)dudv  =  J  J  Gp(u,v)‘3F(t~u,n~v)  dudv, 


at  Rt 


Rr 


U 


30D*F(t,n)  _  0p  *  F  (  t ,  n ) .  The  second  half  of  the  theorem  fol’ows 


at 


at 


from  the  proof  of  Theorem  2.1. 

Remark . 

Theorem  2.2  could  have  been  phrased  as  follows:  If  ^F  is 

at 

3  F 

piecewise  continuous  and  m  <  3F  <  M,  then  m  <  0D  *  —  <  M. 

at  at 

Thus,  at  points  where  2F_  is  not  continuous,  the  convolution  of  F 

at 

by  Gp  "blends"  the  first  partials  of  the  grid  in  one  section  into 
those  of  the  next.  This  blending  also  occurs  for  all  higher 
derivatives . 

Theorem  2.3. 

If  F(t»n)  is  piecewise  continuous,  then  Gp  *  F(t,n)  is  C°°. 
Proof . 

By  change  of  variables  we  have 
3Qp*F(t,n) 


at 


£-=•  II  0p(  u  ,  v)  •  F(  t-u ,  n-v)dudv 
dt  R  ^ 

Rp 

n+p.  t+p 

3 _  J  J  0  (t+u,  n+v)  F(u,v)dudv 

H  n-p  t-p  p 

n+p  t+p 

/  I  3  0p( t+u,  n+v)  ,  ,  , 

n-p  t-p  — _  •  F(u/V)duav. 

at 


i  • 


tw- 


>V.-1 


r 


f  'i 


<-  .**  .  * 


Since  Op  is  (f  ,  Op  *  F  is  C  . 

Theorem  2.4  ensures  that  if  T  is  on6-to-one  and  p  is  "small1 
then  Qp  *  T  is  one-to-one. 

Theorem  2.4. 

If  T(^,n)  =  (X(£»n)»  Y{£,n))  ‘'here  X(C,n)  and  Y  ( (; ,  n  )  are  C* , 
then  if  the  Jacobian  J  of  T  is  not  zero  there  is  a  p>  0  such 
that  the  Jacobian  Jp  of  0p  *  T(£,r,)  =  (  Op  *  X(£,n),  Op  *  Y  ( £ ,  n )  ) 
is  not  zero.  , 

Proof . 


The  Jacobain  Jp  =  (0p  *  3X)  •  ( 0p  *  3Y)  -  (0p  *  3X)  •  (0p  *  3Y)  . 

at  an  an  at 

By  Theorem  2.2  lim  Jp  =  J. 

p  +0 

Theorem  2.5. 

Let  T(  £ ,  n  )  =  (X( t ,n )  ,  Y(t,n)}  be  a  cl  transformation  and 

let  y  be  the  angle  between  u  =  (AX,  AX  )  and  v  =  ( 3X  f  a_Y  ). 

at  at  an  an 

If  Yp  denotes  the  angle  between  up  =  (0p  *  AX  ,  0p  *  BY  )  and 

at  at 

VP  =  (Op  *  AA  ,  Op  *  AX  >  *  then  lim  =  Y  . 

3n  3n  p^O 

Proof . 


By  Theorem  2.2  lim  up  =  u  and  lim  vp  =  v.  Since  cos  fp  =  _Hp 

p-*-0  p->0  |  Up 


Theorem  2.6. 


If  FU,n)  =  FiU/H)  +  F2U,n)  /  where  and  F2  are  piecewise 
continuous,  then  Qp  *  F(£,n)  =  Op  *  F^(£,n)  +  Op  *  F2(£,n). 

Proof . 

The  additivity  of  the  integral  is  all  that  is  needed  to 
prove  this  theorem. 

Theorem  2.1. 

If  F(£»n)  =  f(C)  -g( n)  where  f  and  g  are  piecewise  continuous, 
then  0p  *  F(£,n)  =(  0p  *  f  ( 5 ))  *(  0p  *  gin))* 

Proof . 

This  result  is  an  immediate  consequence  of  Fubini's  Theorem. 

Theorem  2.8. 

if  FU,n)  =  a  +  BE,  +  cn  +  D^n , then  Gp  *  FU,n)  =  F(C,n). 

Proof . 

It  is  a  routine  check  to  show  that  if  G(£,n)  =  £ >  then 
Op  *  G(£,n)  =  G  (  £  ,  n  )  .  The  theorem  is  now  an  immediate  consequence 
of  Theorems  2.6  and  2.7. 

Theorem  2.8  shows  that  if  a  portion  of  the  grid  of  the  air¬ 
craft  can  de  defined  by  equations  of  the  form  F(£,n)  =  A  +  B£,  +  Cn  +  D£n 
hen  the  convolution  step  can  be  bypassed.  Thus,  the  grid  can  be 
generated  much  more  rapidly  in  that  region. 


III.  The  Examples 


Theorem  3.1. 

If  T{  E,,  n)  -  (X(£,  n),  Y(£,  n)  )  has  piecewise  continuous 

partials  which  satisfy  m^  8_X  <  m..?  c[X  <  M..,  m^-^  <  3Y  <  M^, 

8C  3n  3 C 

.  / 

and  m_  <  _9Y  5  M__  for  points  in  (u-p,  u  +  p]  x  [v-p,  v  +  p],  then 


>  0 


the  Jacobian  Jp  of  Op  *  T  is  non-zero  if  either  n'11m22  _  M21M12 


or  M11M22  -  m^ 2m21  <  °* 

Proof . 

By  the  remark  after  Theorem  2.2  m^  <  Op  *  9  X  ( g  ,n  )  <  M2i' 

etc.  Since  m^^m22  ~  M12M21  -  Jp  -  M11M22  ~  m12m21'  ;*'S  DOt  zero' 

Let  -p  =  xo  <  xi  <  • ■ • <  x2n  =  P  a  partition  of  [-p,p] 
such  that  x^ =  h  for  all  i.  Let  yj  =  xj.  Let  A^j  =  {(i,j)|i,j/0 
or  2n  and  either  i  or  j  odd}.  Let  =  {(i,j)|  i  or  j  =  0  or  2n}. 

Theorem  3.2.  (Simpson's  Rule  in  2-Dimensions) 


If  F(u,v)  is  a  C  function  on  RD,  then 


JJ  F ( u , v ) dudv  =  y3h 2 -1  F(xi,  yj)  -  +  2/3h2*[  F(x£,yj)  +  Err, 


where  |  Err  |  £  M-h  /8  •  (2p)  .  (The  constant  M  =  max  {  3_  F  t  9  F 

3  3  u  ^  ■  3  v 1 

94F  }.) 

,.4 


Proof  : 

Except  for  the  fact  that  a  2-dimensional  version  of  Taylor's 
Theorem  is  necessary,  the  proof  of  this  theorem  is  the  same  as  the 
familiar  Simpson's  Rule. 

Since  the  functions  F(u,v)  to  be  integrated  in  this  paper  are 


X(u,v)  or  Y(u,v)  convolved  with  Op  ,  F(u,v)  =  0  whenever  u  =  +  p 

or  v  =  +  p.  Thus,  for  the  purposes  of  this  paper  JJ  F(u,v)dudv 

RP 

4/3h2*£  R(x.,y.).  For  the  applications  illustrated  in  this 
Aij  1  D 


paper  h  was  chosen  to  be  p/4.  This  choice  of  h  is  equivalent  to 


9. 


dividing  the  square  into  16  equal  pieces.  To  approximate  the 
integral  of  F  when  h  =  p/4  the  function  must  be  evaluated  at  40 
points.  See  Figure  5. 

The  following  example  has  been  worked  out  to  illustrate  how 
the  theory  from  the  provious  section  can  be  applied  to  a  specific 
transformation.  The  equations  given  below  approximate  the  projec¬ 
tion  of  one  half  of  the  X-24C  aircraft  into  the  plane.  The  trans¬ 
formation  to  be  smoothed  is  T(u,v)  —  (X(u,v),  Y(u,v))  where  X(u,v) 
and  Y(u,v)  are  defined  below. 

If  u  —  1,  X(u,v)  =  u  and  Y(u,v)  =  2v. 

If  l<u<2,  X(u,v)  =  -13  +  14u  +  .8u  -  .  8uv  and 

Y ( u , v )  =  -3.5  +  3 . 5u  +  1 . 9v  +  (,l)uv. 

If  2  <  u<  3,  X ( u , v )  =  7  +  4u  -  1.15v  +  .175uv  and 
Y ( u , v )  =  - . 5  +  2u  +  1 . 3v  +  . 4uv . 

Let  F ( u , v )  = Jyj u2  +  v2  +  u  . 

Let  A  =  22  -  /2/4F{l,v)3  +  3/2/4 *v *F ( -1  ,v  )  , 

B  =  76  -  19u  -  2 . 5v  +  .625uv, 

C  =  5.5  -  /2  /  4  F  (  - 1 ,  v  )  3  +  3/2/4 -v-F(  1  ,v)  ,  and 
D  =  22  -  5.5u  +  lOv  -  2.5uv. 

If  3 <  u<  4,  X(u,v)  =  ( u-3  5  A  +  B  and  Y(u,v)  =  (u-3)C  +  D. 

If  u>4,  X ( u , v )  =  22  -  /2/4 -F(5-u,v)3  +  3  /T/  4*v-F(u-5,v)  and 

Y(u,v)  =  5.5  -  /2/4F  ( u-5  ,  v  )  3  +  3  *  v  *  F  (  5-u ,  v ) /2/4  . 

Note  that  the  transformation  on  the  region  u  a  4  is  nothing  more 

% 

than  a  translated  and  rotated  version  of  the  map  w  =  z  .  For 
3  -  u  s  4  the  transformation  is  an  interpolation  between  a  straight 

line  and  the  w  =  z  2  map.  The  grid  for  this  transformation  is 

/ 

illustrated  in  Figure  4.  Note  the  singularities  of  the  first 


derivative  along  the  lines  u=l,  u=2,  u  =  3,  and  u=4.  A  smoothed 


version  of  this  grid  is  seen  in  Figure  5.  This  grid  was 
by  convolving  T  with  the  bump  function  Op(u,v)  where  p  = 
Note  that  the  singularities  have  now  disappeared. 


generated 

0.5. 


The  mathematics  in  this  paper  shows  that  it  is  possible  to 
patch  a  grid  together  from  local  grids.  Even  if  the  "patched" 
grid  is  not  smooth,  it  can  be  approximated  by  a  smooth  one. 

Desirable  properties  such  as  orthogonality  and  appropriate 
clustering  of  grid  points  will  be  almost  retained.  While  this 
approximation  technique  will  never  produce  grids  as  "perfect" 
as  those  generated  by  conformal  or  hyperbolic  techniques,  it 
should  be  useful  in  piecing  together  complicated  configurations 
where  one  is  more  interested  in  obtaining  a  "reasonable"  grid 
rather  than  a  flawless  one.  Once  the  equations  for  the  grid  in 
each  section  have  been  obtained,  the  method  is  very  fast.  Also, 
since  the  convolutions  are  evaluated  locally  there  is  no  accumulated 
error,  (i.e.  Errors  incurred  in  calculating  one  grid  point  do  not  enter 
into  the  calculating  of  the  next.)  The  examples  illustrated  were 
all  run  single  precision.  Obviously,  if  more  accuracy  is 
warrented,  the  computer  programs  could  be  run  with  double  precision 
and  with  a  larger  selection  of  points  when  applying  Simpon’s  Rule. 

When  applying  the  smoothing  techniques  indicated  in  this  paper, 
care  must  be  used  when  choosing  the  value  of  the  parameter  p.  If 
p  is  too  large  the  approximation  will  not  be  close  enough.  Thus, 
the  Tacobian  could  become  zero  or  the  control  on  orthogonality 
could  be  lost.  If  p  is  small  relative  to  the  number  of  grid  points, 
the  grid  will  have  "numerical  discontinuities"  in  the  derivatives. 

One  final  remark  should  be  made.  The  parameter  p  does  not 


have  to  be  a  constant.  If  very  tight  fit  of  grid  lines  is  needed 
at  some  point,  the  parameter  p  can  be  allowed  to  approach  zero. 

This  new  convolution  will  still  be  C  if  p  is  constant  near  points 
of  discontinuities  of  the  partial  derivatives  of  the  transformation 
If  p  varies  arbitrarily  the  new  transformation  will  only  be  . 
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SMOOTHING  PATCHED  GRIDS 
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Abstract 

The  purpose  of  this  paper  is  to  indicate  how  smoothing  techniques  can  be 
utilized  in  the  area  of  grid  generation.  The  focus  of  the  paper  is  to  show 
how  one  global  grid  can  be  patched  together  from  a  number  of  smaller  ones. 

The  procedure  usually  takes  place  in  two  steps.  First,  one  global  grid  is 
patched  together  from  a  number  of  smaller  ones,  allowing  for  the  possibility 
that  the  derivatives  along  common  boundaries  may  not  be  continuous.  The 
second  step  is  to  then  approximate  this  grid  by  a  smooth  one  in  such  a  way 
that  the  essential  structure  of  each  patch  is  preserved. 


This  research  was  supported  in  part  by  AFOSR  Grant  083-0158.  The  paper 
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I.  Introduction 

The  author  developed  the  Ideas  outlined  in  this  paper  while  attempting  to 
construct  a  grid  for  an  aircraft.  A  plane  can  be  divided  in  a  natural  way 
into  components  such  as  the  forebody,  the  airfoil,  the  tail,  etc.  The  regions 
surrounding  these  components  can  usually  be  subdivided  in  a  natural  way  so 
that  suitable  local  grids  (or  patches)  can  be  found  for  each  subregion. 

However,  the  frustration  of  trying  to  splice  together  these  various  pieces 
into  one  global  grid  drove  the  author  to  search  for  a  technique  which  would 
blend  one  patch  into  the  next  while  still  preserving  the  essential  structure 
of  each  local  grid.  The  principal  smoothing  technique  described  here  involves 
a  convolution  of  the  grid  transformation  with  a  "bump"  function  to  obtain  a 
new  smoother  grid  which  approximates  the  old  one.  Each  new  grid  point  can  be 
thought  of  as  a  weighted  average  of  nearby  points. 

In  this  paper  no  effort  will  be  made  to  deal  with  patches  that  overlap  as 
Steger,  Dougherty,  and  Benek  [1]  have  done.  In  fact  the  standing  assumption 
will  be  that  adjacent  patches  will  have  common  boundary.  Moreover,  the  grid 
points  on  a  common  boundary  between  two  adjacent  regions  will  be  assumed  to 
agree.  In  the  terminology  of  M.  M.  Rai  [2]  the  grid  may  have  metric  discon¬ 
tinuities  but  no  discontinuities.  Figure  1  indicates  the  difference.  (The 
author  would  prefer  to  say  the  grid  is  continous  but  not  smooth.) 

Actually,  the  problem  of  smoothing  grids  has  been  encountered  before. 

For  example,  the  elliptic  method  [3]  or  [4]  can  be  thought  of  as  a  smoothing 
technique.  The  reason  for  this  is  that  before  the  iterative  scheme  is  to 
begin,  the  user  must  provide  an  Initial  guess  (smooth  or  not) .  With  good 
forture  this  guess  is  then  rapidly  molded  into  a  smooth  grid.  Even  simpler  is 
the  Laplace  operator 

Xn(I,J)  -  [Xn_L (I ,  J  -  1)  +  X^d,  J  +  1)  +  Xn_x(I  -  1,  J)  +  X^d  +  1,  J)]/4 
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A  few  iterations  with  this  operator  and  a  grid  can  be  smoothed  signif¬ 
icantly.  However,  too  many  iterations  may  lead  to  a  grid  which  is  not 
one-to-one.  Examples  illustrating  this  difficulty  are  discussed  in 
Section  III.  Kowalski  [5]  has  developed  a  variation  of  the  Laplace  operator 
to  smooth  an  algebraic  grid.  He  allowed  his  operator  to  sweep  through  the 
grid  as  many  as  12  times. 


II .  The  Smoothing  Theory 

In  order  to  explain  the  examples  presented  in  Section  III  it  is  first 
necessary  to  present  the  background  to  the  smoothing  theory.  While  grid 
generation  is  primarily  concerned  with  a  discrete  set  of  lattice  points,  it 
will  be  convenient  here  to  present  the  theory  in  terms  of  continuous  func¬ 
tions,  derivatives,  and  integrals.  The  transition  from  the  continuous  theory 
to  the  discrete  theory  will  be  explained  in  Section  III. 

To  develop  the  smoothing  theory  it  is  first  necessary  to  explain  the  term 
"bump"  function.  If  R  denotes  the  real  numbers  and  p  >  0,  then  a  nonnega¬ 
tive  function  0  :  R  +  (0,»)  is  a  bump  function  supported  on  the  interval 
I  — p , p]  if  it  is  smooth,  is  identically  zero  outside  [~p,p],  and  has  the 
property  that 


P  /*" 

0  (x)  d x  «  I  0  (x)  dx  *  1. 

p  p  J -<*  P 

While  there  are  innumerable  choices  for  a  bump  function,  the  one  indicated 

below  is  convenient  to  explain  and  use. 

First  define  the  function  a:R  [0,®)  by  the  rule 


a(x) 


x  <  0 
x  >  0 


Note  that  a  is  C°*  and  identically  zero  on  (-®,0]  .  The  graph  of  a  is 
indicated  in  Figure  2.  If  Bp(x)  -  a(x  -  p)  •  a(p  -  x) ,  then  0p  is 
nonnegative,  is  identically  zero  off  the  interval  [~p,p] ,  and  is  C“.  The 


graph  of  6p  is  indicated  in  Figure  3. 


rp 

If  A  *  8(x)  dx  and 

P  J- P  • 


0  (x)  -  8  (x)/A  ,  then  0  (x)  is  a  Lump  function  supported  on  [ — p , p] . 


Since  grid  generation  is  primarily  concerned  with  arrays  in  2  and  3 

dimensions,  the  notion  of  bump  function  must  be  extended  to  the  square 

R  ■  [-p,p]  x  t“P»Pl*  Note  that  6  (x,y)  ■  0  (x)  •  0  (y)  is  nonnegative, 
P  P  P  P 


identically  zero  off 


R  ,  and  is  C°*. 
P 


Note  also  that 


0p(x,y)  dx  dy 


is 


-  1. 


o 

Definition:  If  F(£,n)  is  a  piecewise  continuous  function  from  R  *  R,  then 
the  convolution  of  F  by  dp  is  defined  by  the  equation 


6p  *  FU.n) 


0  (u,v)  du  dv  . 
P 


Intuitively,  the  convolution  can  be  thought  of  as  an  average  of  the 
values  of  F  over  the  square  [£  -  p,  £  +  p]  x  [n  -  p,  n  +  p]  relative  to 
the  weight  function  0p  (u,v).  In  particular,  if  F(x,y)  ■  1  for  :all 
(x,y)  e  R2,  then  0p  *  F(£,r\)  *  1  for  all  £,n*  Thus,  it  becomes  clear  that 
if  F  is  nearly  constant  near  (£,n),  then  0p  *  F(£,n)  is  "close”  to 
F(£,n)«  Theorem  2.1  gives  a  precise  statement  of  this  observation.  In  fact 
Theorems  2. 1-2. 9  give  precise  formulations  of  the  following  statements. 

1.  If  p  is  "small”,  then  0p  *  F  is  "close"  to  F. 

2.  For  any  p,  the  convolution  0p  *  F  is  a  C"  function. 

3.  If  F(£,n)  is  differentiable  at  (£,n)  and  p  is  small,  then  all  the 
derivatives  of  0p  *  F(£,tO  will  be  close  to  the  derivatives  of 
F(C,n). 


4. 

If 

T(£,n) 

“  (X(£,n),  Y(£,n)),  the  Jacobian 

J 

of 

T  is  not  zero. 

and 

p  is 

small,  then  the  Jacobian  Jp  of 

9P 

*  T 

is  not  zero. 

5. 

If 

T  is 

orthogonal  and  p  is  small,  then 

0P 

*  T 

is  "nearly" 

orthogonal. 

6.  The  convolution  operator  is  linear. 
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7.  The  convolution  operator  is  invariant  when  applied  to  functions  of 


the  form  F(x,y)  ■  A  +  Bx  +  Cy  +  Dxy. 

At  this  point  the  reader  who  is  not  interested  in  the  theory  can  skip 
ahead  to  the  examples  in  Section  111.  An  excellent  reference  for  the  Advanced 
Calculus  needed  in  the  proofs  of  the  following  theorems  is  Rudin  [6] . 

Theorems  2. 1-2. 3  are  lifted  directly  from  the  Differential  Topology  books  by 
Hirsch  [7]  and  Munkres  [8]. 


Theorem  2.1  If  F  is  continuous,  e  >  0,  and 

|F(e  -  u,  n  *  v)  -  F(£,n)]  <  c  for  all  (u,v)  e  Rp,  then 

I  6p  *  F(£, n)  -  F(£,n)  |  <  £  • 

Proof 

The  result  follows  from  the  following  sequence  of  equalities  and 
inequalities. 


|e  *  F(c,n)  -  FU,n)| 


I  /  /  F(5  -  u,  n  -  v)  0p(u,v)  dudv  -  F(e,n)| 


'//« 

R 

P 

R 


v)  •  |f( C  -  u,  n  “  V)  -  F(C,n)|  dudv 


<  e 


//vu- 

R 


v)  dudv  ■  e« 


Theorem  2.2  If  F(£,n)  is  piecewise  continuous,  then  0  *  F(C,q)  is  C*. 


T.  r  i  •  , 


Proof 


By  change  of  variables  we  have 


3e  *  F(5,n) 


~  JJ  ep(u,v)  .  F(5  -  u,  T)  -  v)  dudv 


a  re+p 

—I  /  0 (C  +  u,  n  +  V)  F(u,v)  dudv 

35  Jn-P  -Vp  P 

rn+p  /-E+p  30  (5  +  Uf  n  +  V) 

I  I  - E - - -  •  F(u,v)  dudv 

•'n-n 


0  (£  +  u,  n  +  v)  F(u,v)  dudv 
P 


rrp  5“P 


Since  0p  is  C  ,  6p  *  F  is  C  . 

Theorem  2.3  shows  that  if  p  in  "small"  and  '^If2^  is  continuous  at 

^0  *  F(5,n)  3FfF  v  5 

(C,n)»  then  — — — -  is  "close"  to  — - Note  that  Theorem  2.3  can  be 

easily  generalized  to  higher  derivatives. 

.  30_  *  F(C,n)  aFf.  v 

Theorem  2.3  If  F(C,n)  is  C  ,  then  —2 - — -  ■  9  *  — - 

-  at  p 

Moreover,  if  e  >  0  and  [ _u>  ^ <  e  for  (u,v)  e  R„, 

30  *  F(?,n)  3F(C,h) 

then  1-^-^ - 3^-1  <  c- 


Proof 


Since  JJ 0p(u,v)  •  F(£  -  u,  n  -  v)dudv  -  JJ 8p(u,v)  • 


3F(C  -  u,  n  -  v) 


dudv. 


30„  *  F(c,n) 


■  9  The  second  half  of  the  theorem  follows  from  the 

P  3€ 


proof  of  Theorem  2.1. 


Remark 


Theorem  2.3  could  have  been  phrased  as  follows:  If  is  piecewise 

u,  3P 

continous  and  m  <  —  <  M,  then  m  <  9  *  -r-r  <  M.  Thus,  at  points  where 

P 


t-  '.-  v 1 j.  v  .•  ■-•Vv'-j'- -  .'-.'.'x.  .^V-V>Y*V-V- 
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—  is  not  continuous,  the  convolution  of  F  by  0  "blends’*  the  first 
3£  P 

partials  of  the  grid  in  one  section  into  those  of  the  next.  This  blending 
also  occurs  for  all  higher  derivatives.  Theorem  2.4  ensures  that  if  T  is 
one-to-one  and  p  is  “small",  then  0p  *  T  is  one-to-one. 


Theorem  2.4  If  T(£,n)  -  (X(5,n),  Y(£,n))  where  X(£,n)  and  Y(£,n)  are 
C1 ,  then  if  the  Jacobian  J  of  T  is  not  zero  there  is  a  p  >  0  such  that 
the  Jacobian  Jp  of  0p  *  T(£,n)  -  (0p  *  X(£,n),  6p  *  Y(C,n))  is  not  zero. 


Proof 


The  Jacobian  J 

P 

By  Theorem  2.3. 


Theorem  2.5  Let  T(£,n)  *  (X(£,n),  Y(£,n))  be  a  C1  transformation  and  let 
V  be  the  angle  between  u  *  -||j,  and  v  -  ,  -|^j.  If  Fp  denotes 

the  angle  between  up  -  (@p  *  -g  ,  0p  *  g)  and  v?  «  (op  *  gj  ,  0p  *  g). 

then  lim  T  ■  T. 

p-*0  P 


Proof 

By  Theorem  2.3  11m  u  ■  u  and  lim  v 

p-*0  p  p-*0  P 

and  cos  f  -  *-  T-V  ,  lim  f  -  ?. 

Iul  *  M  p-*0  p 


v.  Since  cos  T 
P 


Theorem  2.6.  If  F(£»h)  *  F^(£,n)  +  F2(£»n)»  where  F^  and  ?2  are 
piecewise  continuous,  then  0p  *  F(£,ti)  ■  0p  *  F^(C,n)  +  0p  *  F2(5*n)» 


Proof 

The  additivity  of  the  integral  is  all  that  is  needed  to  prove  this 


theorem. 


Theorem  2.7  If  F(£,n)  -  f(£)  •  g(n),  where  f  and  g  are  piecewise  contin¬ 
uous,  then  0  *  F(£,n)  “  [0„  *  f(5)l  •  [fl_  *  g(n)]  • 

D  1  p  P 


Proof 


This  result  is  an  immediate  consequence  of  Fubini' s  Theorem. 


Theorem 


2.8  If  F(£,n)  -  A  +  B?  +  Cn  +  D£n,  then  0p  *  F(£,n)  -  F(£,n) 


Proof 

It  is  a  routine  check  to  show  that  if  G(£,n)  ■  5,  then 
0p  *  G(5,n)  *  G(£,n)*  The  theorem  is  now  an  immediate  consequence  of 
Theorems  2.6  and  2.7. 

Theorem  2.8  shows  that  if  a  portion  of  the  grid  of  the  aircraft  can  be 
defined  by  equations  of  the  form  F(£,n)  ■  A  +  B£  +  Cn  +  D£n,  then  the  convo¬ 
lution  step  can  be  bypassed.  Thus,  the  grid  can  be  generated  much  more 
rapidly  in  that  region. 

The  next  theorem  is  of  interest  because  it  gives  sufficient  conditions 
which  ensure  that  the  Jacobian  will  be  nonzero  at  (S,n)  even  in  the  case 
that  the  partial  derivatives  of  T(£,n)  are  not  defined  at  (g,n).  In  most 
applications  these  hypotheses  will  be  satisfied. 


Theorem  2.9  If  T(£,n)  ■  (X(£,n)*  Y(5,n))  ^as  piecewise  continuous  partlals 

and  ad  -  be  >  0,  where  a  <  ~,  <  b,  <  c,  and  d  <  ■—  for  all  points 

35  3h  35  oh 

in  (u  -  p,  u  +  p)  x  (v  -  p*  v  +  p) ,  then  the  Jacobian  Jp  of  0p  *  T  is 


nonzero. 


By  the  remark  after  Theorem  2.3  a  <  0  *  — ,  0  *  ■£—  <  b,  0  *  -r^  <  c, 

P  3£  P  3n  P  3C 


av 

and  d  <  0  *  —  so  that 

P  3h 


-  (eP  *  -S)  -  (eP  *  -  (%  *  *  (S  *  -H)  >  -  -  -  >  »■ 


One  final  remark  should  be  made  at  this  point.  While  all  the  theorems  in 
this  section  have  been  stated  in  a  two-dimensional  setting,  each  one  general¬ 
izes  to  three  dimensions. 


III.  The  Examples 


The  purpose  of  this  section  Is  to  Indicate  how  the  theory  from  Section  II 
can  be  applied  to  smooth  a  patched  grid.  Immediately ,  we  are  confronted  with 
four  problems. 

1.  The  size  of  the  parameter  p  must  be  fixed. 

2.  The  bump  function  must  be  selected. 

3.  A  numerical  integration  technique  must  be  chosen. 

4.  A  method  must  be  found  to  fix  the  points  on  the  boundary,  while 
retaining  the  overall  smoothness  of  the  rest  of  the  grid. 

In  the  discrete  setting  the  point  (£,n)  must  be  replaced  by  the  lattice 
point  (I,J),  where  I  and  J  are  integers.  Since  the  choice  of  the 
parameter  p  determines  the  size  of  the  square  Rp,  the  selection  of  p  now 
becomes  a  decision  concerning  the  number  of  neighbors  of  (I,J)  to  be  used  when 
convolving  with  the  bump  function.  If  the  point  (I,J)  is  to  be  in  the  center 
of  the  square,  the  reasonable  choices  seem  to  be  9,25,49,  etc.  In  the  dis¬ 
crete  situation  the  larger  the  number  of  points,  the  smoother  the  new  grid 
will  be.  However,  the  increased  number  of  computations  could  easily  become 
prohibitive.  For  the  purposes  of  the  examples  presented  here  the  author  chose 
49  neighbors  for  each  point  (I,J).  These  49  will  be  of  the  form 
(I  +  K,  J  +  K')  ,  where  |K|  <  3  and  |K' j<  3. 

While  the  choice  of  bump  function  is  important,  it  does  not  seem  to  be  as 
critical  as  some  of  the  other  problems.  However,  if  one  is  careless  and 
chooses  a  bump  function  which  is  very  near  zero  except  at  the  point  (I,J), 
then  very  little  or  no  smoothing  will  take  place.  Since  the  bump  function  is 
identically  zero  on  the  boundary  of  the  square  Rp,  the  integration  is  now  to 
be  performed  over  the  square  (I  +  K,  J  +  K'),  where  |K|  <  2  and  (K* |  <  2. 
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The  method  of  integration  chosen  is  the  two-dimensional  version  of  the 
Newton-Cotes  formula  indicated  in  Proposition  3.1.  This  two-dimensional  inte¬ 
grator  was  obtained  by  discarding  the  remainder  term  and  taking  the  tensor 
product  with  itself. 


Proposition  3.1  (See  page  93  of  Hildebrand  [9].)  If  xq  <  x^  < 


<  V 


h  *  x  j  -  xt,  and  f  -  ^x^),  then 


I 


6  f(x)  dx  -  yjg-  (41fo  +  216fx  +  27 f 2  +  272f3  +  27f4 


+  216f5  +  41£6>  '  TIW  f<8)(2)- 


In  the  theory  outlined  in  Section  II  there  is  no  mention  of  boundary 
points.  If  the  problem  is  ignored,  the  surface  of  the  aircraft  will  be 
smoothed  along  with  the  grid.  Sharp  corners  will  become  rounded.  (Compare 
Figs.  4  and  5.)  While  there  are  a  variety  of  ways  to  deal  with  this  problem, 
the  author  chose  to  reduce  the  amount  of  smoothing  for  grid  points  near  the 
boundary  by  linear  blending.  In  particular,  if  T(I,J)  denotes  the  original 
grid  point  and  d  is  the  distance  from  T(I,J)  to  the  boundary,  then  the  new 
grid  point  will  be  Tn(I,J)  *  dTs(I,J)  +  (1  -  d)  T(I,J),  where  TS(I,J) 
denotes  the  corresponding  point  of  the  smoothed  grid.  Figures  6-10  indicate 
the  rough  and  smooth  versions  of  various  shapes.  Since  the  value  of  d 
was  never  quite  equal  to  1  in  these  examples,  small  discontinuities 
in  the  derivatives  were  propagated  into  the  interior  of  the  regions. 
Despite  this  problem  the  grids  were  still  smoothed  fairly  well.  In 
the  future  the  author  plans  to  develop  a  distance  function  which  is 
exactly  equal  to  1  throughout  most  of  the  interior  of  the  region  so 
t-hat  the  grid  will  be  smooth  away  from  the  boundary. 
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Figures  11  and  12  have  been  included  to  compare  the  Laplace  operator  and 
the  simple  average  of  the  nine  immediate  neighbors.  While  these  two  methods 
are  both  much  faster  than  the  convolution,  they  usually  need  to  be  iterated  to 
be  effective.  When  iterated  without  boundary  control,  the  grid  points  may 

drift  outside  the  region  in  question.  This  phenomenon  is  demonstrated  by 
Figure  13.  This  same  Figure  also  shows  that  a  transformation  which  is  a 

solution  of  Laplace’s  equation  may  fail  to  be  one-to-one  even  if  it  is 
one-to-one  on  the  boundary.  If  boundary  control  is  forced  on  the  Laplace 
operator,  then  the  grid  points  will  stay  in  the  region.  However,  discon¬ 
tinuities  in  the  derivatives  may  appear  near  the  boundary  as  indicated  in 
Figure  14.  Figure  15  illustrates  the  result  when  the  convolution  operator  is 
iterated  twice.  Obviously,  any  further  iterations  and  the  grid  will  become 
overlapping. 


IV.  Concluding  Remarks 


The  mathematics  In  this  paper  shows  that  It  Is  possible  to  patch  a  grid 
together  from  local  grids.  Even  If  the  patched  grid  is  not  smooth,  it  can  be 
approximated  by  a  smooth  one.  Desirable  properties  such  as  orthogonality  and 
clustering  of  grid  points  will  be  almost  retained.  These  techniques  can  be 
thought  of  as  postprocessors  to  remove  discontinuities  in  the  grid  derivatives 
along  the  boundaries  of  adjacent  patches.  A  future  application  could  be  the 
smoothing  of  grids  created  from  patches,  where  one  patch  is  generated  by  a 
hyperbolic  method,  a  second  by  an  elliptic  method,  a  third  by  an  algebraic 
method  (see  Ref.  [10]  or  [11]),  etc.  The  final  grid  would  then  be  a  smoothed 
version  of  Che  union  of  the  patches. 

One  further  remark  seems  to  be  in  order.  When  the  author  began  this 
research,  he  only  considered  the  convolution  operator.  At  the  suggestion  of 
P.  Eiseman  the  Laplace  operator  was  also  considered.  While  both  worked  well, 
the  Laplace  operator  is  easier  to  program,  faster  in  terms  of  CPU  time,  and 
seemed  to  generate  somewhat  smoother  grids.  A  reason  for  using  the  convolu¬ 
tion  operator  is  that  it  seems  to  be  better  at  preventing  overlapping  near  the 
boundary.  The  Laplace  operator  can  frequently  be  iterated  successfully 
5-10  times.  The  convolution  works  best  when  applied  1-2  times. 
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